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Two linear least squares test problems, both fifth decree polynomials, have been run on more than 
twenty different computer programs in order to assess their numerical accuracy. Among the programs 
tested were representatives from various statistical packages as well as some from the SHARE library. 
Essentially five different algorithms were used in the various programs to obtain the coefficients of the 
least squares fits. The tests were run on several different computers, in double precision as well as 
single precision. By comparing the coefficients reported, it was found that those programs using 
orthogonal Householder transformations or Gram-Schmidt orthonormalization were much more accu- 
rate than those using elimination algorithms. Programs using orthogonal polynomials (suitable only 
for polynomial fits) also proved to be superior to those using elimination algorithms. One program, 
using congruential methods and integer arithmetic, obtained exact solutions. In a number of programs, 
the coefficients reported in one test problem were sometimes completely erroneous, containing not even 
one correct significant digit. 
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1 . Introduction 

Since the time when the electronic computer began to supplant the desk calculator as the 
chief tool for solving linear least squares problems, numerous least squares computer programs 
have been written. These programs have utilized a variety of computational algorithms. Because 
least squares problems are by their very nature frequently ill-conditioned, the numerical accuracy 
achieved by a least squares program strongly depends upon the choice of the algorithm. Many 
programs have been written which use methods appropriate for desk calculators but inappropriate 
for computers. Anscombe II] 1 has aptly remarked: "Textbooks of statistical method display a won- 
derful unanimity in recommending computational procedures that are suited to desk calculators 
but are perilous for computers. Only with some determination can the statistician break himself of 
bad habits and become adequately informed about round-off error." 

The present study was undertaken to assess the numerical accuracy of representative least 
squares programs from a variety of sources. Two test problems, both fifth degree polynomials, 
have been run on more than twenty different programs. Included in the study were programs from 
the BMD Biomedical Computer Programs collection, the C-E-I-R Multi-Access Computing Services 
library, the IBM SHARE library, the IBM System/360 Scientific Subroutine Package, the Univac 
MATH-PACK and STAT-PACK collections, and the Project MAC 7094 disk files. A detailed list- 
ing of the sources of the programs is given in appendix A, together with a brief description of each 
program. 

For a number of programs, the test problems were run in double precision as well as in single 
precision. This, of course, necessitated certain changes in the original programs. 

The programs included in this study used essentially five different algorithms: orthogonal 



1 Figures in brackets indicate the literature references at the end of this paper. 
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Householder transformations, Gram-Schmidt orthonormalization, orthogonal polynomials, Gauss- 
ian or Jordan elimination, and a congruential method with computations in integer arithmetic. 
Previous studies appraising linear least squares programs and comparing the results of dif- 
ferent algorithms have been made by Cameron [9], Freund [20], Bright and Dawkins [7],Zellner and 
Thornber [46J, Longley [29], and Jordan [27]. The present study differs from the earlier ones mainly 
by including a larger selection of widely used and easily accessible programs. 

The linear least squares problem may be briefly stated as follows: One has n observations or 
measurements of a "dependent" variable y which are statistically independent with common vari- 
ance cr 2 whose expected values are given by a linear function of the corresponding values of k 
"independent" variables, xi 9 #2, . . ., Xk, k ^ n. In matrix notation we say that the n observations 
have expected values E(Y) =Xfi, where Y is an n X 1 vector, X is an n X k matrix, and ft is a k X 1 
vector of unknown coefficients. Assuming that X is of rank k, the least squares estimates of the 
coefficients are given by J3 = (X f X)~ l X'Y. Other quantities of interest are Y = Xf3, the vector of pre- 
dicted values; 8=Y—¥, the vector of residuals; and s 2 = (Y— f)'(Y— Y) , an estimate of 

1 n — k 

the variance cr 2 . 

In running certain programs, modifications were occasionally made to input and output 
formats. Other changes were made in five of the programs using elimination algorithms because 
the original versions of these programs failed to give solutions to the fifth degree polynomial prob- 
lems. In particular, features that may have been intended to prevent execution of computations 
subject to excessive rounding error were sometimes bypassed. Details of these changes, and some 
remarks on the effectiveness of the features which were bypassed, will be given in section 8. 

Four computers were used: the GE 235, the IBM 7094, and the Univac 1107 and 1108. The 
1108 which was used is located at the National Bureau of Standards, and the 7094 which was chiefly 
used is located at Harry Diamond Laboratories, Washington, D. C. The programs run on the 235, 
the 1107 and the Project MAC 7094 utilized consoles at the National Bureau of Standards connected 
to computers at other locations. 



2. The Test Problems 

The two main test problems which were used throughout this investigation are identified as 
Y\ and Yl. Both were fifth degree polynomials, with the values of x being the integers 0, 1, 2, . . ., 
20. The "observations," Y\ and F2, were calculated from the following equations: 

Yl: y=l + x + x 2 + x 3 + x 4 + x\x = 0(l)20, 

Y2: y=l + 0.lA:-h0.01x 2 + O.OOU 3 -h0.0001x 4 + 0.00001^ 5 , * = 0(1)20. 

Thus the values of Yl were integers having from one to seven digits, and those of Y2 were five- 
decimal numbers ranging from 1.00000 to 63.00000. 

If the least squares solutions were computed with no rounding error, one would obtain 



0(11)- 



/3(Y2) = 



1. 
0.1 

,01 

.001 

.0001 
. .00001J 



and for both problems the residual standard deviation would be zero. 

For some programs the input required was the 21 values of x andy. Some programs required, 
in addition, the powers x 2 , x 3 , x 4 , and x 5 to be entered as input. Other programs required as input 
the 6 by 6 matrix X'X and the 6 by 1 vector X'Y. It should be noted that the elements olX'X are 
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integers having from 2 to 14 digits, the elements of X'Y for Fl are 8- to 14-digit integers, and the 
elements of X'Y for F2 are 5-decimal numbers having up to 13 significant digits. The input is listed 
in table 9. 

The two test problems, Fl and F2, were chosen because they are so highly ill-conditioned that 
some programs fail to obtain correct solutions while other programs succeed in obtaining reason- 
ably accurate solutions. Polynomial problems were chosen because polynomial fitting is an im- 
portant type of linear least squares problem which occurs frequently in practice. 

The ill-conditioning of the two test problems can be described more explicitly. One measure 
of the condition of a matrix A is the P-condition, defined as 



/ 



P(A) = 



where A. is the numerically largest eigenvalue of A and \x is the numerically smallest eigenvalue of 
A. (See Newman [34, p. 240]). 

For A=X'X, the 6X6 matrix associated with Fl and F2, the P-condition is 4.095 X 10 13 . In 
this respect, it is similar to the Hilbert matrix of order 10, whose P-condition is L603 X 10 13 (see 
Fettis and Caslin f IT]). The P-condition of the Hilbert matrix of order 11 is 5.231 X 10 14 . The relation 
between the Hilbert matrix and the matrix X'X which arises in a polynomial fit is discussed in 
Forsythe [18]. 

Most of the programs which were tested obtained more accurate solutions for F2 than for Fl. 
If we let A denote the 7X7 matrix 



A = 



X'X X'Y 
Y'X 



we find that for F2, P(A) =4.095 X 10 13 , whereas for Fl, P(A) =6.829 X 10 13 , indicating that the 
system involving Fl is more ill-conditioned than that involving F2. 

The test problem used by Longley [29] was also highly ill-conditioned. For the 7X7 matrix 
X'X of his problem, the P-condition is 2.361 X 10 1! ' 

3. Summary of the Results 

Tables 1 to 6 present a brief summary of the main results. A count, C ; , of the number of correct 
significant digits in each computed coefficient was obtained as follows: 

Let fij (j= 1,2,. . ., 6) denote the "true" value of the coefficient — that is, the value computed 
with no rounding error. Let /3j denote the value calculated by the computer. Then 



ogio 



Pj-Pj 



if Ift-jSjI^Oandfr^O 



Cj = { - tofto|ft-&|, if |ft-&| * and (3j = 

Z), the approximate number of decimal digits with which the machine computes, if Pj — J3j = 0. 

The above approach to counting the number of correct digits in a computed value has been 
used by Jordan [27] and others. 

Tables 1 to 6, in the columns headed "Average Number of Correct Digits" report 

1 6 

From the above definition, a negative count can occur. For example, if j8/= 1.0, and J3j= 136.0, 
we get Cj = — 2.130. This indicates that J3j is wrong by roughly two orders of magnitude. 
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Table 1. Summary of programs run in single precision — 8 digits 



Program 


Com- 
puter 


Algo- 
rithm a 


Average number of 
correct digits 


Rank 


Date 


of run 


Fl 


F2 


Fl 


F2 


Fl 


F2 


ALSQ 


1108 
1108 
7094 
1108 
7094 
1108 
7094 
1108 
1108 
7094 
7094 
1108 
7094 
1108 
1108 
1108 
1108 
1108 
1108 
1108 
1108 
7094 


HT 
E 
E 
E 
E 
E 
? 

HT 
OP 

E 
E 
E 

GS 

GS 

GS 

GS 

E 

E 

E 

E 

E 

E 


4.098 

-0.106 

0.742 

-0.123 

1.389 

-0.264 

-2.756 

4.528 

2.118 

-0.140 

-0.607 

-0.907 

3.954 

4.137 

-1.976 

3.593 

-0.191 

-0.658 

0.066 

0.066 

0.651 

-5.300 


5.368 
1.981 
1.721 
2.287 
2.312 
2.622 

-0.301 
5.840 
4.363 
1.856 
1.460 
1.224 
5.968 
5.464 
0.419 
6.197 
2.280 
1.527 
2.767 
2.767 
2.407 

-2.871 


3 
12 

8 
13 

7 

16 
21 

1 

6 
14 
17 
19 

4 

2 
20 

5 
15 
18 
lC* 
lOi 

9 
22 


5 

14 
16 
12 
11 

9 
21 

3 

6 
15 
18 
19 

2 

4 
20 

1 
13 
17 

7i 

71 
10 
22 


10-18-67 

12-13-67 

12-30-66 

12-18-67 

4-12-67 

3- 5-68 

5-17-68 

5- 1-68 

4-12-68 

5-16-67 

12- 9-66 

2-29-68 

12- 5-66 

10-18-67 

3- 7-68 

9-10-68 

10- 7-68 
11-14-67 

11- 7-67 
11- 8-67 

7- 2-68 
6-28-67 


10-18-67 


BMD02R 


11-17-67 


BMD03R... 


1- 3-67 


BMD03R... 


12-18-67 


DAM 


4-12-67 


DAM 


3- 5-68 


UNFIT (Miller). 


8-15-68 


LSTSQ : 


5- 6-68 


MATH-PACK ORTHLS... 


4-12-68 


MPR3 


5-18-67 


OMNITAB (Invert) 


12- 9-66 


OMNITAB (Invert) 


2-29-68 


OMNITAB (Ortho) 


12- 5-66 


OMNITAB (Ortho) 


10-18-67 


ORTHO (no iteration) 


3- 7-68 


ORTHOL 


9-10-68 


POLRG 


10- 7-68 


SPVMTX 


11-14-67 


STAT-PACK, GLH 


11- 7-67 


STAT-PACK, REBSOM 


11- 9-67 


STAT-PACK, RESTEM 


7- 2-68 


WRAP 


6-28-67 







a E = Elimination method; 
OP = Orthogonal polynomials. 



GS = Gram-Schmidt orthonormalization; HT = Orthogonal Householder transformations; 



TABLE 2. Summary of programs run in single precision — 9 digits 



Program 


Com- 
puter 


Algo- 
rithm a 


Average number of 
correct digits 


Rank 


Date of run 


Fl 


F2 


Fl 


F2 


Fl 


F2 


LINFIT*** 


235 
235 
235 
235 
235 
235 
235 


E 
E 
GS 
OP 
E 
E 
E 


0.905 
0.308 
4.102 
3.349 
1.402 
0.612 
1.169 


2.894 
2.483 
6.354 
5.922 
3.213 
2.920 
3.183 


5 

7 
1 
2 
3 
6 
4 


6 

7 
1 
2 
3 
5 
4 


12- 1-67 

12-28-66 

1-25-67 

2-19-68 

12-30-66 

11-30-67 

1- 3-67 


12- 1-67 


LSCF— *** 


12-28-66 


LSFITW*** 


1-25-67 


POLFIT 


2-19-68 


SIMEX-*** 


1- 5-67 


STAT20*** 


11-30-67 


STAT21** 51 : 


1- 3-67 







a E = Elimination method; GS = Gram : Schmidt orthonormalization; OP = Orthogonal polynomials. 



Table 3. Summary of programs run in double precision— 16 digits 



Program 


Com- 
puter 


Algo- 
rithm a 


Average number of 
correct digits 


Rank 


Date of run 




Fl 


F2 


Fl 


F2 


Fl 


F2 


BMD05R 


7094 
1107 


E 

E 


6.953 
7.882 


6.230 
9.959 


2 

1 


2 
1 


1- 5-67 
1-23-67 


1- 5-67 


DPVMTX 


1-23-67 







a E = Elimination method. 
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Table 4. Summary of programs run in double precision — 18 digits 



Program 



Com- 
puter 



Algo- 
rithm 



Average number of 
correct digits 



71 



F2 



Rank 



Yl 



Y2 



Date of run 



Yl 



Y2 



ALSQ 

BMD02R 

BMD05R 

DPVMTX 

LSTSQ 

MATH-PACK, ORTHLS 

ORTHO 

ORTHCHno iteration).... 

ORTHOL 

POLRG 

STAT-PACK, RESTEM. 



1108 
1108 
1108 
1108 
1108 
1108 
1108 
1108 
1108 
1108 
1108 



HT 

E 

E 

E 

HT 

OP 

GS 

GS 

GS 

E 

E 



12.667 

9.645 

9.368 

9.744 

14.643 

12.098 

13.188 

7.963 

13.212 

9.290 

9.494 



15.322 
12.865 
11.791 
13.484 
16.293 
14.461 
15.514 
10.354 
15.604 
11.806 
12.019 



4 
7 
9 
6 
1 
5 
3 

11 
2 

10 
8 



10-19-67 
4-17-68 
9-10-68 
2-27-68 
7-22-68 

10-16-68 
1-29-68 
3- 7-68 
9-25-68 

10- 7-68 
7- 1-68 



1-29-68 
4-17-68 
9-10-68 
2-27-68 
7-22-68 

10-16-68 
1-29-68 
3- 7-68 
9-25-68 

10- 7-68 
7- 1-68 



a E= Elimination method; GS= Gram-Schmidt orthonormalization; HT= Orthogonal Householder transformations; 
OP= Orthogonal polynomials. 

Table 5. Summary of programs run in single precision (8 digits) with inner products accumulated in double precision 

(18 digits) 



Program 


Com- 
puter 


Algo- 
rithm ;| 


Average number of 
correct digits 


Rank 


Date of run 




Y\ 


K2 


Y\ 


Yl 


Y\ 


Yl 


ALSQ 


1108 
1108 
1108 


HT 
HT 

GS 


3.506 
8.000 
3.904 


6.530 
6.279 
6.459 


3 
1 

2 


1 
3 
2 


10-18-68 

4-30-68 

10-21-68 


10-18-68 


LSTSQ 


5- 6-68 


ORTHO... 


10-21-68 







a GS = Gram-Schmidt orthonormalization; HT= Orthogonal Householder transformations. 



Fable 6. Summary of program run in multiple precision integer arithmetic 





Program 


( !om- 

puter 


Algo- 
rithm a 


Rational form 
Floated form 


Average number 
of correct digits 


Date o 


f run 




Yl 


Y2 


Yl 


Y2 


SOLVER 


1108 


C 


X 

18.000 


oc 
17.347 


7-16-68 


7-16-68 







a C = Congruential method. 

For two programs reported in table 1, BMD03R run on the 7094 and DAM run on the 7094, 
the count for several coefficients was made in a different manner. The BMD03R program printed 
the coefficients in a fixed-decimal format, with five decimals. The DAM program used a floating- 
point format with only three decimals printed. A coefficient printed as 0.00010, when the true 
coefficient was 0.0001, was given a count of 2, and 0.100E01, when the true coefficient was 1., 
was given a count of 3. In such cases the assigned count may have been too small, since the coeffi- 
cients may have been calculated accurately to more digits than were printed. In running these 
two programs on the 1108, the output format was changed so that eight significant digits were 
printed. 
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Each of the tables (1 through 6) summarizes a set of results for a particular machine preci- 
sion— 8, 9, 16, 18, etc. digits. Within each table the various programs are ranked for each of the 
two test problems, with rank 1 denoting the best performance according to the count C. 

Table 1 includes single precision (8-digit) programs run on two different computers, the 7094 
and the 1108. It was felt that combining the results from two computers was justified in view of the 
similar performance of the four programs which were run on both computers. These four programs 
were BMD03R, DAM, OMNITAB (using INVERT), and OMNITAB (using ORTHO). The average 
number of correct digits for the two test problems from the 1108 agrees with the corresponding 
average from the 7094 to within 0.9 digits except in the case of Fl run on DAM, where the difference 
is 1.653. This larger difference may possibly be attributable to modifications in the program. Further- 
more, other test problems have been run on OMNITAB (using ORTHO) on both computers, and 
again the results were quite similar. 

The symbols in the Algorithm column of the tables denote the following: 
C Congruential method, integer arithmetic 
E Elimination method 

GS Gram-Schmidt orthonormalization 

HT Orthogonal Householder transformations 

OP Orthogonal polynomials. 

From time to time at a given computer installation changes are made in hardware and software 
with the result that a particular job run on two different days may not produce identical numerical 
output. For this reason the date of the computer runs is included in the tables. 

Details (individual coefficients and counts) supporting tables 1 through 6 are given in ap- 
pendix B. 

4. Programs Using Orthogonal Householder Transformations 

LSTSQ is a program written by Peter A. Businger using orthogonal Householder transforma- 
tions. This algorithm is described by Golub [21] and Businger and Golub [8J. The program applies 
a sequence of orthogonal transformations to the n by k least squares matrix X to obtain a decomposi- 
tion X = QR, where R is upper triangular and Q'Q = h-. A pivoting strategy is used so that at each 
step the column with the largest sum of squares is reduced next. Once an initial solution is obtained, 
the program iterates to obtain a (possibly) improved solution. 

Of all the programs using floating-point arithmetic included in this study, LSTSQ appears to 
have given the best performance. In table 4, which reports the performance of eleven double preci- 
sion programs, we see that LSTSQ ranked first for both Fl and F2. In table 1, which reports the 
performance of 22 single precision programs, we see that LSTSQ obtained rank 1 for Fl and rank 
3 for F2. Ranks 1 and 2 for the F2 problem were obtained by ORTHOL and OMNITAB (using 
ORTHO), two programs using Gram-Schmidt orthonormalization which will be discussed more 
fully in the next section. Table 5 reports the performance of three programs which used single preci- 
sion arithmetic except for the accumulation of inner products, where double precision arithmetic 
was used. Here we see that LSTSQ ranked first for Fl (having a perfect score of 8.000) and ranked 
third for F2. In the two instances just mentioned where LSTSQ ranked third for F2, we see that 
the difference separating it from the top-ranking program is small, being 0.357 in table 1 and 0.251 
in table 5. 

Golub and Businger recommend that all inner products be accumulated in double precision. 
By comparing tables 5 and 1 we see that when LSTSQ included this feature, the average counts 
increased from 4.528 to 8.000 for Fl and from 5.840 to 6.279 for F2. With all operations performed 
in double precision (see table 4), the counts increased to 14.643 and 16.293, respectively. 

The other program using Householder transformations was ALSQ, written by G. W. Stewart, 
III. This program contains no pivoting and no iteration. In tables 1, 4, and 5 we see that ALSQ 
performed not quite as well as LSTSQ which included these features, except in one instance. In 
this one instance, F2 in table 5, we note that its performance was slightly better than that of LSTSQ. 
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By examining tables 1 and 5, one may see the effect of accumulating inner products in double 
precision versus accumulating them in single precision. As one would expect. ALSQ did better 
in computing the coefficients of F2 when the double precision accumulation was included. For 
Fl, surprisingly, we see that ALSQ lost accuracy with this feature included. A look at the details 
of the program revealed how this phenomenon occurred. After the matrix X has been decomposed 
to obtain QR (as described earlier), the coefficients are computed by back-substitution. The first 
coefficient to be computed is /3e (the coefficient for the fifth-degree term), and this is obtained 
from one arithmetic operation, a division. Correctly calculated to ten digits, this division is 
21011.77901 



21011.77901 



= 1. In the single precision version, the coefficient was calculated as 

& = 2101l'713 = L000000 ° (t<) 8 si g nificant di S its ) (1) 

whereas in the version with inner products accumulated in double precision the calculation was 

P" = 21011 753 = LOOO W04 (to 8 significant (,i ^ its )- (2) 

We note that in (1), both the numerator and denominator in question are farther from their true 
values than in (2), but they are closer to each other, so that j§<; in (1) happens to be closer to the 
true value of /3<;. Subsequently, /3<; enters into the calculations of the five other coefficients with 
the result that all the coefficients for Fl from the single precision version are slightly more ac- 
curate than those from the version using double precision inner products. 

5. Programs Using Gram-Schmidt Orthonormalization 

ORTHO is a program written by Philip J. Walsh using a Gram-Schmidt orthonormalization 
process. This algorithm is described by Davis and Rabinowitz [13|, [14|, Davis [12], and Walsh 
[42]. ORTHO exists as a FORTRAN program, an ALGOL procedure, a BASIC program, and as a 
routine of the OMNITAB program (see Hilsenrath, et al., [23]), where it is called by the commands 
FIT and POLYFIT. 

Starting with the n X k matrix X, the Gram-Schmidt process of ORTHO obtains ^p — XT ,x 
and/3 =T'~ l <p'Y, where T' ~ x is upper triangular and <p'(p = h. This algorithm includes a feature of 
reorthonormalizing the vectors of (/?. proceeding from a first approximation <f>j to a (usually) better 
approximation <pj. From table 1 it is clear that this reorthonormalizing is vital to the algorithm, for 
ORTHO's good performance in handling Fl and Y2 deteriorated when this iteration was omitted. 
For Fl, the count of correct digits dropped from 4.137 to— 1.976. and for Y2 the drop was from 5.464 
to 0.419. In table 4, also, we see that in double precision the omission of the iteration resulted in 
a loss of about five correct digits for both problems. 

The LSFITW*** program, written in BASIC, listed in table 2 also uses the ORTHO algorithm. 
The computer used for running the programs of table 2 works with about one more decimal digit 
than the computers covered in table 1, so one would expect more accuracy from LSFITW*** 
than from OMNITAB (ORTHO). We find slight improvement for F2, but no improvement for Fl. 
Of the seven programs reported in table 2. LSFITW*** ranked 1 on both problems. Note that 
there are no Householder transformation programs included in table 2. 

The ORTHO program was also run in a version using single precision except for the accumu- 
lation of inner products, where double precision was used. In table 5 we see that there were three 
programs in this category, and for both problems ORTHO ranked second. Its performance on 
F2 improved by about one digit compared to the performance of the ORTHO version using strictly 
single precision. On the Fl problem, however, there was a slight loss in accuracy. Actually, three 
coefficients gained accuracy and three lost accuracy, with a net loss in the average count. (A 
similar loss which occurred with ALSQ was discussed in the previous section.) In ORTHO, the 
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final calculation to obtain the coefficients /3 is the matrix multiplication j8= (T')~ l a, where (7") _1 
is an upper triangular matrix such that (T')- 1 T- 1 = (A'A) -1 , and a = T l X'Y. Nearly all of the 
nonzero elements of (T')~ 1 and a were more accurately computed when inner products were 
accumulated in double precision than when this feature was omitted. In the three coefficients of 
Y\ which lost accuracy, an examination of the details showed that in the individual multiplica- 
tions involved in the matrix multiplication, the version using double precision for inner products 
was always more accurate than the strictly single precision version. But in the final addition of 
the various products, where the terms have alternating signs, there was heavy cancellation and the 
errors combined in such a way that the jSj's from the single precision version happened to be closer 
to the true values of the coefficients than those computed with double precision accumulation of 
inner products. 

ORTHOL is a program using a modification of the Davis-Rabinowitz algorithm written by 
James W. Longley and Roger A. Blau [30], It differs from Walsh's ORTHO in two respects: 
(1) the iteration procedure includes the dependent variable as well as the independent variables, 
and (2) before any other operations are applied to the matrix X, from each element of each vector 
of X. the truncated mean of that vector is subtracted. (The "truncated mean" denotes the largest 
integer less than or equal to the mean if the mean is nonnegative, and the smallest integer greater 
than or equal to the mean if the mean is negative.) ORTHOL obtained the top rank for F2 in single 
precision, but ranked fifth for Y\ (table 1). In double precision (table 4), it ranked second on both 
problems. 

6. Programs Using Orthogonal Polynomials 

Since the two test problems are both polynomial fits, we were able to test programs in which 
the algorithm used orthogonal polynomials. This method, described by Forsythe [18], is attractive 
because it generally requires many fewer operations than other methods. 

Two such programs were included in this study. One was the UNIVAC 1108 MATH-PACK 
routine, ORTHLS (see Programmers Reference [40]). The other was POLFIT, an anonymous 
program written in BASIC. 

In tables 1,2. and 4 we see that the performance of the orthogonal polynomial programs is 
not as good as that of the Householder transformation and the Gram-Schmidt programs (with itera- 
tion), but the performance is better than that of any of the programs using elimination algorithms. 
This finding is in agreement with the results of Bright and Dawkins [7] who ran a number of poly- 
nomial test problems via two methods: matrix inversion using a Gauss-Jordan reduction, and 
orthogonal polynomials. In all cases they found the orthogonal polynomial method superior. 

7. A Multiple Precision Integer Arithmetic Program Using Congruential Methods 

Morris Newman, in his paper "Solving Equations Exactly" [35] described a congruential 
method for finding the exact solution of a system of linear equations Ax = b where the elements of 
A and b are all integers. His FORTRAN program SOLVER will solve systems in which A is a square 
matrix at most 100 by 100 and the elements of A and b are numerically less than 10 20 . This method 
is not at all sensitive to the condition of A, but it can be time-consuming for large systems. 2 The 

solution is printed in two versions: -(1) jc= ( -j j) z, where z is a vector of integers and the deter- 
minant det A is an integer, and (2) x in floated double-precision format, accurate to about 17 digits 
on the 1108 

The two test problems, Fl and Y2, were run on this program, as indicated in table 6. The 
input required was the matrices XX and X'Y. Since the elements of X'Y for F2 are not integers, 
it was necessary to multiply these numbers by 100,000 before obtaining the solution. 



2 The running time on the Univac 1108 for the solution of VT and F2 was 11 seconds, including 5 seconds for compilation of the program. The six problems de- 
scribed in the latter part of this section, all having 6x6 systems, required 14 seconds, including 4 seconds for compilation. To solve a 20 X 20 system, the worst pos- 
sible case requires about 30 seconds, and an "average" case takes less time. For a 40 X 40 system, an average case requires about one minute, and the worst possible 
case requires about six minutes. A "bad" 100 X 100 case might require 40 minutes or more running time. 
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Having a program which produces exact solutions, we can determine what will happen to the 
solution when we round the input, the elements of X X andX'Y. These elements for Yl are integers 
having no more than 14 digits. Six additional problems were run in which the input was succes- 
sively rounded to 13, 12, 11, 10, 9, and 8 significant digits. The effect of this rounding was to change 
the solution dramatically. In table 7, which gives the solutions to these six problems (rounded to 
10 decimals), we see that "small" changes to the elements of X'X and^T produce "large" changes 
to the solution, ft. 

At first glance, the fact that the coefficients calculated for the problem having input rounded 
to 13 significant digits agree with the coefficients obtained from unrounded input to only 4, 5, 6, 
or 8 digits seems quite surprising. But with a few simple calculations we can see why the agreement 
is no better than it is. First, referring to table 9, we note that there is only one 14-digit number in 
X'X, and since this ends with a zero, rounding to 13 significant digits leaves X'X unaltered. In 
X'Y only the last element has 14 digits. Here, 25,537,373,767,266 was rounded to 25,537,373,767,270. 
Let B= (bij), i, j— 1, . . .,6, denote the inverse oiX'X and let (qj)=X'Y,j= 1, . . .,6. Consider 
)8i . the first coefficient. We have 

<* 5 

0i = b n q\ H- hvzq-i + 613^3 + b u q 4 + ftistfs + &i6<?6 = ^ &ij<# + bieqe- 

j 1 

The only quantity which is affected by changing from unrounded data to rounded data with 13 
significant digits is q^. Now 

6,« = -28,046,715,376,452,025,326,796,800/(Det X'X), where 

De\X'X= 1,677,193,579,511,831,114,542,448,640,000. Since q 6 = 25,537,373,767,266 we have 

6,«g 6 = ~ 716,239,453,512,581, 907,404 ,696,441, 196,473,548,800/(Det X'X). 
Also, 

2 b xj qj= 716,239,455, 189,775 ,486,916,527,555,738,922,188,800/(Det X'X). 

In combining the last two numbers, one positive and the other negative, we see that the first 
eight digits of the two numerators are identical, so that the numerator of /3i has eight fewer digits 

than has V b\jqj or bisqe- 

In solving for fi\ from input rounded to 13 significant digits, we have ^- = 25,537,373,767,270 
so that 

buw=~ 716,239,453 ,512,694,094,266,202,249 ,297,780,736,000/(Det X'X). 

The numerator here differs from the numerator of &i«qr« in the thirteenth significant digit, but 

5 

after combining this with ^ b\jqj and losing eight significant digits, we obtain 

5 
#= X bijqj+bi6q« 

_ 1,677,081,392,650,325,306,441,141,452,800 
1,677,193,579,511,831,114,542,448,640,000 

= 0.99993 31104 (to 10 decimals). 

There are similar losses of significant digits in calculating the other coefficients. 

A rigorous presentation of the sensitivity of the solution of a system of equations Ax = b with 
respect to variations in A and b is given in Wilkinson [43, p. 91] and Wilkinson [44, p. 189]. 
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Table 7. Exact solutions to approximate problems 

The column A(N) gives a count of the number of digits in the solution of (X'X)/3 = X'Y for input rounded to N digits which are in agreement with the solution for un- 
rounded input. The unrounded input (elements of X'X and X'Y) consisted on integers having no more than 14 digits. 



Solution for unrounded 


Solution for input 


A(13) 


Solution for input 


A(12) 


input 


rounded to 13 sig. digits 




rounded to 12 sig. digits 






0.99993 31104 


4.175 


0.99672 54481 


2.485 




1.00015 06443 


3.822 


1.00718 28792 


2.144 




0.99994 21662 


4.238 


0.99727 93936 


2.565 




1.00000 79679 


5.099 


1.00037 12813 


3.430 




0.99999 95470 


6.344 


0.99997 90427 


4.679 




1.00000 00091 


8.043 


1.00000 04168 


6.380 




Avg. = 


= 5.287 


Avg.= 


= 3.614 




Solution for input 


A(ll) 


Solution for input 


A(10) 




rounded to 11 sig. digits 




rounded to 10 sig. digits 






1.06051 32874 


1.218 


0.91988 69708 


1.096 




0.86927 95602 


0.884 


1.19460 48731 


0.711 




1.04911 27057 


1.309 


0.92297 46106 


1.113 




0.99333 63623 


2.176 


1.01080 85953 


1.966 




1.00037 44589 


3.427 


0.99937 78148 


3.206 




0.99999 25798 


5.130 


1.00001 25555 


4.901 




Avg.= 


= 2.357 


Avg.= 


2.166 




Solution for input 


A(9) 


Solution for input 


A(8) 




rounded to 9 sig. digits 




rounded to 8 sig. digits 






3.70810 09327 


- 0.433 


-24.56199 35653 


-1.408 




-4.34362 06781 


-0.728 


52.60541 27451 


-1.713 




2.92257 14823 


- 0.284 


- 17.89853 20445 


- 1.276 




0.74656 29295 


0.596 


3.52648 11096 


- 0.403 




1.01394 46989 


1.856 


0.85941 47174 


0.852 




0.99972 81121 


3.566 


1.00276 62377 


2.558 




Avg.= 


= 0.762 

i 


Avg.= 


-0.232 



8. Programs Using Elimination Algorithms 

The majority of the programs tested in this investigation used some form of an elimination 
algorithm. Although this was the most popular method, it was the least successful. None of these 
programs performed as well as those using Householder's transformations, Gram-Schmidt ortho- 
normalization (with iteration), or orthogonal polynomials. 

Within this class of programs, there were several variations in the method of obtaining the 
least-squares coefficients. In some cases, the matrix XX was inverted, after which the inverse 
was postmultiplied by X'Y to obtain (3= (X'X)- l X'Y. In one program the matrix 



XX 





X'Y 

1 



was inverted. Here, the inverse is 



A~ 



(X'xy 
o 



. Another program inverted the matrix Z'Z 



where the vectors of Z were obtained from the vectors ofX by subtracting the mean of each vector 
from every element of that vector. A number of programs obtained the solution by inverting a matrix 
of correlation coefficients. The five stepwise regression programs made use of matrix partitioning 
in connection with inverting a matrix of correlation coefficients. 

8.1. Stepwise Regression Programs 

The five stepwise regression programs were BMD02R, MPR3, the STAT-PACK program 
RESTEM, WRAP, and STAT20***. They all, to a greater or lesser extent, follow Efroymson's 
algorithm [16]. Tables 1, 2 and 4 give the results of these five programs. 
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The BMD02K program is described in BMD Biomedical Computer Programs [15]. The 
UNIVAC 1108 STAT-PACK Program RESTEM is described in the Programmers Reference 
[41]. The two programs, MPR3 and WRAP, are both from the SHARE library [26]. The former 
was written by M. A. Efroymson and the latter was written by M. D. Fimple. The program 
STAT20*** is included in the C-E-l-R Multi-Access Computer Service library, and a brief writeup 
on how to use the program is given in C-E-I-R's "Library Programs Documentation" [10]. 

In running the two test problems on three of the stepwise programs, namely BMD02R, 
RESTEM and STAT20***, calculations stopped before the solutions were obtained. In all three 
programs, computations stopped in the F2 problem after x 9 x 2 , x 4 , x 5 and a constant term had 
entered the regression equation. In the Fl problem, BMD02R and STAT20*** stopped after 
jc 4 , x 5 and a constant term had entered, and RESTEM stopped after x 5 and a constant term had 
entered. In order to obtain the coefficients for the fifth degree equations, certain features of these 
three programs had to be bypassed. 

In the printout of RESTEM and STAT20*** obtained from the original (unaltered) programs, 
there were no clues to indicate that there was a rounding error problem. The initial runs of the 
BMD02R program, however, printed the messages "ERROR TERMINATION IN SQRT ROU- 
TINE" and SQRT CALLED AT SEQUENCE NUMBER 01032 OF MAIN PROGRAM." These 
messages were produced by the computer system, not by the BMD02R program. They indicated 
the nature of some of the trouble — that the argument of a certain square root function was nega- 
tive. The initial BMD02R runs furnished another clue of computational difficulties. The output 
of this program includes various calculated F- values which arc needed for entering and removing 
variables from the regression. In both Fl and K2, there were one or more /^-values (labeled "F TO 
ENTER") which were negative. 

It was found that the RESTEM and STAT20*** programs had also calculated negative F- 
values, and checks involving F-values had to be bypassed in order to obtain the fifth degree solu- 
tions. Moreover, in the RESTEM program it was necessary to change the value of "minus infinity" 
from — 10 : * 8 to — 10 34 before satisfactory results could be obtained for any least squares problem. 

WRAP, the program with the lowest rankings in table 1, computed coefficients which were 
exceptionally far from the true values. These coefficients are listed below. 



Y\ 


K2 


True j$ 


Computed J8 


True/3 


Computed jg 


I. 


2991622. 
- 6065892. 
2218821. 
-296194.5 
16462.20 
-322.5731 


1. 

0.1 
.01 
.001 
.0001 
.00001 


- 33.84546 
71.54880 
-26.16913 
3.493256 
-0.1936966 
.003812985 



Since WRAP performed so poorly on the two test problems, Fl and F2, some other test 
problems were run in order to verify that the program could handle problems which were not so 
badly conditioned. Let (71(A) and U2(k) be defined as follows: 

Ul(k): y=l+x + x 2 +. . .+jc*, 

U2(k): y= 1 + 0.1 * + 0.01 jc 2 +. . . + 10"***. 



Taking x— 0(1)20, k — 1, 2, 3,4, the y-values were calculated for Ul(k) and U2(k). Using these calcu- 
lated y's as input, it was found that the coefficients for degrees 1, 2, and 3 computed by the WRAP 
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program had some accuracy, but those for degree 4 were computed inaccurately. The results for 
degrees 3 and 4 are given below. 



C/l(3) 


C/2(3) 


True/3 


Computed fi 


True/3 


Computed /3 


1. 
1. 
1. 
1. 


1.223297 
0.8286224 
1.022150 
0.9992682 


1. 

0.1 
.01 
.001 


1.000416 
0.09971083 
.01003707 
.0009987662 




171(4) 


£72(4) 


True/3 


Computed /3 


True/3 


Computed ft 




-731.1589 
852.1714 
-193.1382 
15.94353 
0.6330207 


1. 

0.1 
.01 
.001 
.0001 


0.8673919 
.2555575 
- .02548887 
.003731290 
.00003291620 



8.2. Other Programs Using Elimination Algorithms 

Two other BMD programs were tested. The BMD03R program, Multiple Regression with 
Case Combinations, inverts a matrix of correlation coefficients. BMD05R, Polynomial Regression, 
inverts the matrix Z'Z where the vectors of Z are formed from the vectors ofX by subtracting the 
mean of each vector from every element of that vector. All the crucial operations of BMD05R, 
such as the forming of inner products and matrix inversion, are carried out in double precision. 
The performance of BMD03R and BMD05R is shown in tables 1, 3, and 4. 

DAM is a general-purpose computer program for data processing and multiple regression 
written by Rudolf R. Rhomberg, Lorette Boissonneault, and Leonard Harris, International Mone- 
tary Fund [36]. In running the two test problems on DAM on the 1108, computations stopped after 
a fourth degree polynomial was fitted, and the message "INSUFFICIENT NUMBER OF OBSER- 
VATIONS OR DATA ARE ALL ZEROS, PROGRAM CANNOT COMPUTE EQUATION 5" 
was printed. It was found that a computed variance was zero and that this condition causes the com- 
putations to stop. By bypassing the checks on this computed variance, results for fifth degree 
fits were obtained. On the 7094, however, the fifth degree results were reached without any such 
difficulties. DAM's performance on the two computers is given in table 1. 

Two lines of table 1 report the results from OMNITAB on the 7094 and the 1108 where the 
matrix commands INVERT, MMULT, and MTRANS were used. Here the 21 pairs of (x, y) values 
were read into the computer, the powers of x were generated, the matrices X'X and X'Y were 
obtained via MTRANS and MMULT, the inverse of X'X was obtained via INVERT, and p was 
then obtained via MMULT. The solutions were far less accurate than those obtained from OMNI- 
TAB by using the command POLYFIT which calls on the ORTHO routine. 

The program POLRG is the polynomial regression program of the IBM System/360 Scientific 
Subroutine Package [24], [25]. This program calls four subroutines, GDATA, ORDER, MINV, 
and MULTR,in the course of obtaining the least squares coefficients and other quantities of interest. 
These subroutines perform the following operations: 

(1) GDATA generates the powers of the independent variable, finds means and standard 
deviations, and sets up a correlation matrix. 
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(2) ORDER chooses a dependent variable and a subset of independent variables from a larger 
set of variables. 

(3) MINV inverts the correlation matrix using the "standard Gauss-Jordan method." 

(4) MULTR computes the regression coefficients and related quantities, such as the sum of 
squares attributable to the regression and the sum of squares of deviations from the regression. 

We see from table 1 that the single precision version of POLRG obtained rather low scores 
on both test problems. A double precision version of POLRG was also run, and the performance 
here as reported in table 4 was comparable to other programs using similar elimination algorithms. 

The user of POLRG specifies m, the highest degree polynomial to be fitted, and the program 
automatically reports the results of fitting polynomials of successively increasing degrees, starting 
with the first degree. If there is no reduction in the residual sum of squares between two successive 
degrees of polynomials, the program stops the problem before completing the analysis for the 
highest degree specified. In running both test problems, Y\ and Y2, in single precision the analysis 
stopped after degree four, and in lieu of a fifth degree polynomial fit, the message "NO IMPROVE- 
MENT" was printed. In order to complete the calculations for the fifth degree, the checks on 
"improvement" were bypassed. In the double precision version, fifth degree results were obtained 
without any such alterations. 

The Programmer's Manual for the IBM System/360 Scientific Subroutine Package [25] con- 
tains some warnings regarding the accuracy of computations. The reader is informed that the ac- 
curacy of the computations in many of the routines is highly dependent upon the number of signifi- 
cant digits available for arithmetic operations. It is pointed out that matrix inversion and many of 
the statistical subroutines fall into this category, and that the user may, therefore, wish to use double 
precision versions of these routines. (The programs are so constructed that conversion to double 
precision is an easy matter.) An appendix of the manual classifies the subroutines of this package 
into three categories. These are: (1) subroutines having little or no effect on accuracy, (2) sub- 
routines whose accuracy is dependent on the characteristics of the input data, and (3) subroutines 
in which definite statements on accuracy can be made. Only one of the four subroutines called by 
the POLRG program, namely ORDER, is in the first category. The other three subroutines, GDATA, 
MINV and MULTR, fall in the second category. In connection with this second category we read 
that "it cannot be assumed that the results are accurate simply because subroutine execution is 
completed." 

A more explicit statement is given in connection with the subroutine GDATA. Here there is 
a comment in the program stating that if m, the highest degree polynomial to be fitted, is equal to 
5 or greater, single precision may not be sufficient to give satisfactory results. Since the manual's 
test problem for POLRG specifies m=4> and has 15 data points, one might infer that satisfactory 
results would be obtained for this problem. This is not the case, however. In the solution to this 
problem given on page 410 of the manual, the intercept term for the polynomial regression of degree 
4 is reported to be —5.26735. An accurate calculation shows that this term is actually — 6.04262, 
so that the reported term had no correct significant digits. The four reported regression coefficients 
were correctly computed to only one or two digits. Furthermore, the sum of squares of deviations 
from the regression is reported to be 128.85156, whereas it is actually 17.67310. This error is 
also propagated into the calculation of the mean square, the F value, and the improvement in terms 
of sum of squares. The calculated values of Y were found to be correct to one, two or three sig- 
nificant digits, with the residuals correct to one digit or less. 

In concluding this digression concerning the accuracy of the test problem accompanying a 
particular program of a particular package, we note a remark given in the Programmer's Manual 
under "Purposes and Objectives of the Package": "While this package may provide many of the 
tools necessary to solve the more commonly encountered problems in engineering and science, 
there is no intent to imply that these subroutines represent the current state of the art in statistics 
or numerical analysis." 

The programs SPVMTX and DPVMTX appearing in tables 1, 3, and 4 use single and double 
precision versions, respectively, of the same algorithm. These two programs were adapted by 
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Sally T. Peavy, National Bureau of Standards, from two subroutines in the SHARE library: A. R. 
Sadaka's 7090-F1 3180INV1 Single Precision Matrix Inversion with Selective Pivoting and A. R. 
Sadaka's 7090-F1 3181INV2 Double Precision Matrix Inversion with Selective Pivot. Mrs. Peavy's 
adaptations of these programs included accuracy checks on the computed inverse. A brief descrip- 
tion of these accaracy checks is in order. Let A be the input matrix to be inverted and Z the result 
of the inversion. Let E — l — AZ. Let B=(bij) be an n . X n matrix, and let N(B) be a norm denned 
in any of the following ways: 

/ 2\l/2 

N\(B)=l V bij I (the Euclidean or Frobenius norm) 

* '. j ' 

N»(B)=n max |fey-| 
U j 

N 3 (B)=max y 2 N ■ 

In order to guarantee that Z be a good approximation to A~ l , it is only necessary to have 

N(Z) N(E) „ _, . NIZ) ME) . , r , , . . at at a 

— ztttt: — small. 1 his quantity, — , r/ ^, x , is computed lor each of the three norms /Vi, 2V2, and 

l-N(E) l-N(E) 

N3, and provides upper bounds on the error in the elements of the computed inverse. See Newman 
[34, pp. 227— 230J and Taussky [38, pp. 284-286J for a fuller discussion of the norms described 
above. 

The matrix which was inverted by SPVMTX and DPVMTX in solving the two test problems 
was 



XX X'Y 

1 



wnose inverse is 



(ID" 1 -P 
1 



In the single precision solution, the three bounds for Y\ were, respectively, —160. —940, 
and — 140, and those for Y2 were — 2.2, — 7.0 and — 2.7. If an accurate inverse had been obtained, 
the error bounds would have been small positive numbers. That the bounds were negative is a 
clear warning that the computed inverses, including /3, are not accurate. The double precision 
1108 solution obtained the bounds 0.00048, 0.0075, and 0.00046 for Fl, and 0.000000039, 0.00000064, 
and 0.000000087 for F2. In both problems these bounds are quite conservative. In the solution 
of Yl, the largest error in the elements of (X'X)~ l is 5.5 X 10~ 13 , and the largest error in the /3/s is 
5.5 X 10~ 7 . In the solution of F2, the largest error in the elements of {X'X) ~ l is again 5.5 X 10~ 13 , 
and in theft's is3.3XlO~ 13 . 

The 1108 version of OMNITAB now uses the SPVMTX routine for matrix inversion and prints 
out the smallest error bound. OMNITAB reported the smallest error bound for the inverse of 
XX to be —6, a negative number. This was in agreement with the results from inverting X'X 
via the FORTRAN program SPVMTX where the three error bounds were given as —1.7,-6.0 
and -2.1. 

Each of the two STAT-PACK programs, GLH, General Linear Hypotheses, and REBSOM, 
Back Solution Multiple Regression, has its individual features, but for the two test problems the 
solutions were carried out in the same manner, so that the coefficients obtained from the two pro- 
grams were identical, as is indicated in table 1. Both programs invert XX by calling a matrix 
inversion subroutine called JIM which uses a Gauss-Jordan elimination scheme with maximal 
column pivoting and row scaling. The GLH program has an option whereby the user can enter 
restraints in the case X'X is not of full rank. The REBSOM program has the feature that the user 
can enter an F-value to be used as a criterion for removing variables from the regression after an 
initial solution has been computed. 

An error was found in the REBSOM program in the calculation of the variance of Y. After 
estimating k coefficients (including possibly a constant term) from n observations, the formula 

used for the variance of Y is var Y = — ; — . The denominator of this formula should read 

n — k — 1 

n — k rather than n — k — 1. 
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The BASIC program UNFIT***, available in the C-E-I-R Multi-Access Computer Service, 

X'X X'Y' 
Y'X 2y?J 



in order to obtain /3, inverts the matrix A = 

(X'X)- 1 ^- 



whose inverse, if it exists, is 



§L 



Zfi-Y'Y ^y 2 t-Y'Y 



-P' 



1 



J,y\-Y'Y Zfi-Y'Y J 
When Y=Y, the matrix A is singular. In the two test problems Y— F, and the matrix A, if it were 
formed in the computer without any rounding error, would be singular. But A, for Fl and F2, con- 
tains 14-digit numbers, whereas the GE 235 computer works with approximately nine-digit num- 
bers, so that rounding of the elements of A is inevitable, and the version of A contained in the 
computer is not singular. An "inverse" was obtained, and from this fi was immediately computed. 
Table 2 gives the results. 

A third problem was used to test the program LINFIT***. Here we had 



X 



The last column (row) of A is the sum of the three preceding columns (rows), so that A is clearly 
singular, linlike Fl and F2, this problem is such that the elements of A have fewer than nine 
digits. An "inverse" was obtained by LINFIT***, however, and was printed as 



a a a — a 

a a a — a 

a a a — a 

(i — a — a a 



1 5~ 




~6 




~4 


4 12 


20 


1 4 


, Y= 


5 


, and A = 


4 


8 6 


18 


1 2 1 




4 




12 


6 46 


64 


J 2 2 




5 




20 


18 64 


102 



^rhere a =3.67091 X 10< 



To obtain the inverse of A in this program, the matrix command MAT Q=INV(^) is used. 
This command inverts the matrix A by using an elimination method with row pivoting (Kurtz [28|). 
The method is described in section 1.2 of Stiefel [37 j. 

LSCF--*** is another BASIC least squares program available in the C-E-I-R Multi-Access 
Computer Service. Here the solution, /3. is obtained by inverting XX and then post-multiplying 
the inverse by X'Y. Table 2 shows that LSCF--*** ranked below the other programs of this table. 
The inverse of X'X is obtained by using the matrix command INV, the same command as was used 
in LINFIT***. 

The SIMEX-*** program originated at the Naval Ordnance Laboratory. The input required 
for this program was X'X and X'Y , since SIMEX-*** solves n equations in n unknowns. An elimina- 
tion algorithm is used to obtain the solution. The input for this program was limited to nine signifi- 
cant digits. Recalling the results of table 7 which gave the exact solution for the Fl problem when 
input data was rounded to nine digits, it is not surprising that the average number of correct digits 
for Fl, reported in table 2, is only 1.402. This lies between the "accuracy" achieved by the nine-digit 
and ten-digit problems of table 7. 

The BASIC program STAT21*** obtains (X'X)- 1 and )8 by applying Jordan elimination to 
XX and X'Y; the results appear in table 2. 

The LINFIT program included in table 1 is one of eighteen statistical routines described in 
On-line Analysis for Social Scientists by James R. Miller [32]. This library of routines exists in the 
Project MAC 7094 disk files. The two test problems were run on the LINFIT program on a time- 
shared computer via a remote console communicating with Project MAC. A description of Project 
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MAC may be found in Crisman [11]. Miller states that "these routines may be used without exten- 
sive prior training in mathematics, statistics, or computer operations," but in view of LINFIT's 
poor performance on these two problems, it appears that there may be pitfalls in using this program. 
The method used by the L1NFIT program is unknown. By conjecture, it has been included in this 
section among programs using elimination algorithms. 

9. Results from a Problem Having a Nonzero Standard Deviation 

In the two test problems, Fl and F2, treated thus far, the residual standard deviation was 
zero. A third test problem, H*, is one where the standard deviation is nonzero. This problem 
was run on five programs in both single and double precision to see whether the fact that a least 
squares fit has a standard deviation of zero might be a factor influencing the accuracy of 
computations. 

The values of Fl* were derived from the values of Fl by adding 2.0 to the Fl value when x 
is even and subtracting 2.0 from the Fl value when x is odd. The input for Fl* is listed in table 9. 
A fifth degree polynomial fit for the Fl* problem has the solution (to 16 decimals) 



j6(Fl*)- 



2.0459627329192547 
0.1815856777493606 
1.1701440301521066 
0.9870776685960425 
1.0003230582850989 
1. 0000000000000000 



with the residual standard deviation equal to 2.3251684. 

Table 8. Comparison of results from two problems: one with nonzero standard deviation (Yl*) and one with zero standard 

deviation (Yl) 

All problems were run on the 1108 computer 

Single Precision (8 Digits) 



Program 



Algorithm a 



Average number of 
correct digits 



YV 



Yl 



Rank 



n* Fl 



Date of run 



n* 



Yl 



BMD02R 

LSTSQ.... 

MATH-PACK, ORTHLS 

OMNITAB (Ortho) 

POLRG 



E 

HT 
OP 

GS 
E 



0.464 
3.485 
2.053 
3.711 
-0.074 



-0.106 
4.528 
2.118 
4.137 

-0.191 



12-19-67 
10-15-68 
10-15-68 
10-15-68 
10-15-68 



12-13-67 

5- 1-68 

4-12-68 

10-18-67 

10- 7-68 





Double Precision (18 Digits) 












Program 




Algorithm 3 


Average number of 
correct digits 


Rank 


Date of run 




Fl* 


Fl 


Fl* 


Fl 


Fl* 


Fl 


BMD02R 


E 

HT 
OP 

GS 
E 


9.657 

13.913 

12.079 

13.136 

9.270 


9.645 
14.643 
12.098 
13.188 

9.290 


4 
1 
3 
2 
5 


4 
1 
3 

2 

5 


4-17-68 
10-16-68 
10-16-68 

3-27-68 
10-17-68 


4-17-68 


LSTSQ 


7-22-68 


MATH-PACK, ORTHLS 


10-16-68 


ORTHO 


1-29-68 


POLRG 


10- 7-68 







a E= Elimination method; GS= Gram-Schmidt orthonormalization; HT= Orthogonal Householder transformations; OP= Orthogonal polynomials. 
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Table 9. Input for fifth degree polynomials 



X 


Y\ 


Y2 


y 1* 


0. 


1. 


1.00000 


3. 


1. 


6. 


1.11111 


4. 


2. 


63. 


1.24992 


65. 


3. 


364. 


1.42753 


362. 


4. 


1365. 


1.65984 


1367. 


5. 


3906. 


1.96875 


3904. 


6. 


9331. 


2.38336 


9333. 


7. 


19608. 


2.94117 


19606. 


8. 


37449. 


3.68928 


37451. 


9. 


66430. 


4.68559 


66428. 


10. 


111111. 


6.00000 


111113. 


11. 


177156. 


7.71561 


177154. 


12. 


271453. 


9.92992 


271455. 


13. 


402234. 


12.75603 


402232. 


14. 


579195. 


16.32384 


579197. 


15. 


813616. 


20.78125 


813614. 


16. 


1118481. 


26.29536 


1118483. 


17. 


1508598. 


33.05367 


1508596. 


18. 


2000719. 


41.26528 


2000721. 


19. 


2613660. 


51.16209 


2613658. 


20. 


3368421. 


63.00000 


3368423. 



Matrix XX 



21. 

210. 

2870. 

44100. 

722666. 

12333300. 


210. 

2870. 

44100. 

722666. 

12333300. 

216455810. 


2870. 

44100. 

722666. 

1233330O. 

21 64558 10. 

3877286700. 


44100. 

722666. 

12333300. 

216455810. 

3877286700. 

70540730666. 


722666. 

12333300. 

216455810. 

3877286700. 

70540730666. 

1299155279940. 


12333300. 

216455810. 

3877286700. 

70540730666. 

L299155279940. 

24163571680850. 


Matrix AT for Fl 


Matrix X'Y for F2 


Matrix \> for H* 


13103167. 

229558956. 

4106845446. 

74647573242. 

1373802809082. 

25537373767266. 


310.39960 

5058.55410 

87258.40800 

1549291.38666 

28043466.66600 

514843723.46850 


13103169. 

229558976. 

4106845866. 

74647581842. 

1373802985062. 

25537377366266. 



Table 8 summarizes the results, comparing the accuracy of the coefficients for Fl * with 
the corresponding accuracy for Fl. We see that the results for the two problems are quite similar, 
in both single and double precision. The largest differences occurred with the program LSTSQ 
in single precision, where the average number of correct digits was 4.528 for Fl , and the average 
for Fl* was 3.485, a decrease of 1.043. 

On the basis of this comparison it appears that the fact that the standard deviation was zero 
in the test problems did not appreciably affect the accuracy of computations. 



10. Concluding Remarks 



(1) Computational procedures appropriate for desk calculators may be perilous for computers. 

(2) Of the four procedures using floating-point arithmetic which were included in this study, 
orthogonal Householder transformations and Gram-Schmidt orthonormalization proved to ho the 
best. Orthogonal polynomials ranked next. Elimination methods were the least successful but 
the most popular. The multiple precision integer arithmetic procedure using congruent ial methods 
was unique in obtaining exact solutions. 



338-397 0-69— 2 
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(3) Some other algorithms apparently of high quality which have been published in the last 
few years were not included in this study. These include: 

(a) Bauer [2], 

(b) Bjorck and Golub [6] , 

(c) Bjorck [5j. 

Bauer [2] gives an ALGOL procedure using iterative refinement for finding the least squares 
solution of X(3 = K, where X is n X A (A ^ n) of rank A and Y is n X p. The procedure is based on 
the decomposition of X into UDR where U is n X k with orthogonal columns, D = (U'U)~ 1 , and 
R is upper triangular. This decomposition yields a triangular system R/3 = U'Y which is solved by 
back substitution. The reduction to Rf3=U'Y is carried out by a Gaussian elimination scheme, 
but with a suitably weighted combination of rows used for elimination instead of a single row. 

Bjorck and Golub [6] and Bjorck [5] (see also Bjorck [3], [4]) give two least squares algorithms 
with certain common features. Both take advantage of the fact thatZ'8= 0, where 8 is the vector 
of residuals, to obtain the solution /3 in Xfi = Y from the augmented system of n+ k equations: 



/ X 

X' 



M-CT- 



Both algorithms include 8 as well as /3 in the iterative refinement procedure. 

The two algorithms are based on (i) orthogonal Householder transformations, and(ii) a modified 
Gram-Schmidt orthogonalization process. 

Both the classical Gram-Schmidt orthogonalization process and the modified Gram-Schmidt 
orthogonalization process, as described by Bjorck [3], decompose the matrix X into QR where Q'Q 
is diagonal and R is upper triangular. In the classical procedure, at the ith stage, the tth column 
vector is made orthogonal to each of the i — \ previously orthogonalized column vectors; this is 

done for column indices i = 2, 3 A. In the modified procedure which Bjorck uses, at the /th 

stage, the (k — t-hl) column vectors indexed i, i+1, . . . , k are made orthogonal to the '(i — l)-th 
column vector; this is done for column indices i = 2, 3, . . ., A. Jordan [27] shows why the modified 
procedure is superior to the classical procedure. Bjorck [3] states that his modified Gram-Schmidt 
procedure is equivalent to Bauer's method using weighted row combinations mentioned above. 
The algorithms for both the orthogonal Householder transformation method and the modified 
Gram-Schmidt method are generalized to handle the case where X is of less than full rank. In this 
case linear constraints are entered. 

The papers of Bjorck [3, 5] and Bjorck and Golub [6] discuss the number of operations and the 
storage requirements of their least squares algorithms. 

(4) Programmers who have been writing least squares programs, especially for statistical 
packages, have often not been taking advantage of the advances in this area made by numerical 
analysts in recent years. 

(5) The importance of accumulating inner products in double precision cannot be overstressed. 
A number of recent papers on least squares computations have emphasized this point. These 
include Businger and Golub [8], Bauer [2], Golub and Wilkinson [22], Bjorck and Golub [6], and 
Bjorck [5]. On many third-generation computers which have double precision built into the hard- 
ware, double precision arithmetic is quite efficient. 

(6) Iterative refinement is another valuable feature of recent algorithms. The three algorithms 
described in remark (3) above all include this feature. Four programs included in the present 
study (LSFITW***, LSTSQ, ORTHO, and ORTHOL) made effective use of iterative refinement. 
Golub and Wilkinson [22], who discuss this topic, also mention that the condition number of X'X 
is approximately the square of the condition number of X, so that it is advantageous to work with 
X rather than X'X whenever possible. Moler [33] and Forsythe [19] discuss the details of iterative 
refinement in connection with solving nX n systems of linear equations. 
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(7) The users of least squares programs ean take certain precautionary steps to gain an 
awareness of whether or not a rounding error problem exists. Various suggestions were made in 
the previous studies of Cameron, Freund, Zellner and Thornber, and Longley. These 1 suggest ions 
included the following: 

(a) Run test problems where the coefficients are known. 

(b) Transform the data (e.g., by subtracting means). 

(c) Do the calculations several times, scaled differently each time. 

(d) Shuffle the columns of X and run the problem more than once. 

(e) Check whether X'd = 0. 

(f) Use double precision arithmetic. 

(8) Another check on the accuracy of least squares coefficients, suggested by Joseph M. 
Cameron, is the following. After carrying out the usual fit of Y to A independent variables, do a 
second fit, taking F (the predicted values) and refitting to the A original independent variables. 
If there were no rounding error at all, one would obtain exactly the same coefficients from the refit 
as from the original fit, and the standard deviation of the refit would be zero. The extent to which 
the second set of coefficients agrees with the original set can give one some information about the 
severity of rounding error. 

A number of test problems were run on the 7094 and the 1108 in order to investigate the rela- 
tionships among the coefficients of the original fit, those of the refit, and those one would obtain 
if there were no rounding error. The test problems consisted of 55 polynomials with various ranges 
of x, various degrees from 1 to 8, and various coefficients. All 55 were run in both single and double 
precision on the 7094, and twelve of them were run in single and double precision on the 1108. 

In these test problems, the following result was obtained: If the coefficients from the refit 
(denoted by b ) agreed with the coefficients from the original fit (/3) to an average of more than three 
digits, and if the elements of X'X and X'Y can be represented in the computer without rounding, 
then the number of digits in /3 ; in agreement with hj{j= 1, . . ., k) was approximately the same 
as the number of correct digits in 0j. More precisely, whenever the two conditions just stated were 
met, it was found (with one exception) that the number of digits in (3j in agreement with those of 
bj (i)\'m double precision was within 1.0 of the number of correct digits in /3,, and (ii) in single pre- 
cision was within 2.0 of the number of correct digits in /3 h for all j. The one exception occurred in 
a sixth degree polynomial with ,v = — 1()( 1 ) 10. In the double precision run the two sets of coefficients 
agreed to an average of about 12.5 digits, and the elements of X'X and X'Y had at most 13 sig- 
nificant digits; here 0> agreed with S> to 16 digits but was correct only to about 13 digits. 

(9) Efforts to provide comparative data on the amount of computer time required by the 
various programs included in this investigation, as well as comparisons of storage requirements, 
were unsuccessful. The programs which were included in this study originated from many sources, 
and they exhibited considerable variation with respect to what quantities were calculated as well 
as with respect to the methods of calculation. The program ALSQ, for example, at one end of the 
spectrum, calculated simply the coefficients, the residuals, the predicted values, and the residual 
sum of squares for the requested fifth degree polynomial. The single precision version of ALSQ 
required eight seconds to process both test problems on the 1108; the storage requirements were 
709 memory cells for the code and 323 for the data. The double precision version of ALSQ processed 
both problems in seven seconds, requiring 715 memory cells for the code and 618 for the data. 
Nearer the other end of the spectrum was the Biomed program BMD05R. The output here con- 
sisted of the coefficients and their standard deviations for polynomial fits of degrees 1, 2, 3, 4, 5, 
an analysis of variance for degrees 1, 2, 3, 4, 5, predicted values and residuals for degree 5, a 
plot of observed and predicted values for degree 5, and means and correlation coefficients of the 
input data. This program (computing some operations in double precision) required 20 seconds 
on the 1108 to process the two test problems; the storage requirements were 3,119 memory cells 
for the code and 15,168 for the data. It becomes evident that an intercomparison of running time 
among the different programs is not meaningful. 
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Moreover, in repeated runs of a particular program there is fluctuation from run to run in the 
amount of time required. For example, on the same day three separate jobs were submitted to be 
run on OMNITAB (using the command POLYFIT) on the 1108, with the following I results: 

(a) Y\ alone: 8 seconds. 

(b) Y2 alone: 12 seconds. 

(c) Y\ and Y2 together: 8 seconds. 

The 1108 version of OMNITAB requires about 50,000 memory cells for storage. The run times given 
here include unknown components of time for operation of the computer system. 

Although one would expect a double precision version of a particular program to require more 
time than a single precision version, there were several instances on the 1108 where double pre- 
cision required less time than single precision. 

It was outside the scope of this investigation to make a detailed comparison of algorithms with 
respect to efficiency of computation time and storage requirements. Similarly, no comparative 
examination of the outputs provided by the programs was made. Rather, this study focused atten- 
tion on the performance of existing programs. 

(10) In any mathematical calculation carried out on a computer, it is desirable to know whether 
an accurate solution has been obtained or whether the result of a calculation is contaminated by 
rounding error to such an extent that it is worthless. This goal has been achieved in some areas. 
Martin, Peters, and Wilkinson [31], in their paper giving an algorithm for solving Ax — Z>, where 
A is an n X n positive definite matrix and b is an n X p matrix, state that their procedure "either 
produces the correctly rounded solutions of the equation Ax=b or indicates that ,4 is too ill-condi- 
tioned for this to be achieved without working to higher precision (or is possibly singular)." Sim- 
ilarly, Wilkinson's program [45] for the solution of an ill-conditioned system of equations Ax = b, 
where A is rcXn, "gives either a solution of the system which is correct to working accuracy 
or alternatively indicates that the system is too ill-conditioned to be solved without working to higher 
precision and may even be singular." 

It appears that the goal set out above has now been achieved in the linear least squares pro- 
gram of Bjorck and Golub [6]. The authors state that their procedure may be used to compute 
accurate solutions and residuals to linear least squares problems, but that the procedure will fail 
when X modified by rounding errors has less than full rank, and that it will also fail if X is so ill- 
conditioned that there is no perceptible improvement in the iterative refinement. The user is easily 
informed of these situations. 



I would like to express my appreciation to Joseph M. Cameron who suggested this investiga- 
tion and made many valuable contributions, to Joan R. Rosenblatt for helpful discussions, and to 
Joseph Hilsenrath, Russell A. Kirsch, and Thomas Hoover for use of their time-shared computer 
facilities to run several problems. Thanks are due to Gene H. Golub, Stanford University, for his 
constructive remarks. The assistance of James Doyle, Univac Division, Sperry Rand Corporation, 
in debugging one program is also appreciated. 

1 1 . Appendix A. Sources of the Programs, With Brief Descriptions 

ALSQ. A FORTRAN IV subroutine to solve the linear least squares problem, written by G. W. 
Stewart III, Union Carbide Corp., Oak Ridge, Tennessee (present address: The University of Texas, 
Austin, Texas). This program uses a modification of the Businger-Golub algorithm [8]. 
BMD02R, Stepwise Regression. One of the Biomedical Computer Programs, written in FOR- 
TRAN [15]. 

BMD03R, Multiple Regression with Case Combinations. One of the Biomedical Computer 
Programs, written in FORTRAN [15]. 

BMD05R, Polynomial Regression. One of the Biomedical Computer Programs, written in FOR- 
TRAN [15]. 
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DAM. A general purpose computer program for data processing and multiple regression, written 
in FORTRAN by Rudolf R. Rhomberg, Lorette Boissonneault, and Leonard Harris, International 
Monetary Fund [36]. 

DPVMTX. A double precision FORTRAN IV program for inverting a matrix or solving a set of 
linear equations. To a program from the SHARE library (7090-F1 3181INV2 Double Precision 
Matrix Inversion with Selective Pivot, written by A. R. Sadaka [26]), Sally T. Peavy, National 
Bureau of Standards, incorporated accuracy checks. 

LINFIT. A program which fits a linear function to collected data via least squares. Optional 
constraints may be applied to the fitting coefficients to make them nonnegative, add to a constant, 
etc. One of eighteen statistical routines written by James R. Miller [32]. This library of routines 
exists in the Project MAC 7094 in the disk files of user number T169 2750. 

UNFIT***. A program written in BASIC for linear least squares curve fitting and computing 
correlations. Origin: Dartmouth College, Hanover, N.H. Available in the C-E-I-R Multi-Access 
Computer Services library [10]. 

LSCF--***. A least squares polynomial curve fitting subroutine written in BASIC. Origin: Dart- 
mouth College, Hanover, N.H. Available in the C-E-I-R Multi-Access Computer Services library 
[101 

LSFITW***. A least squares curve fitting program written in BASK]. Adapted by John B. Shu- 
maker, National Bureau of Standards, from Philip J. Walsh's ORTHO algorithm [42]. Available 
in the C-E-I-R Multi-Access Computer Services library [10]. 

LSTSQ. A FORTRAN IV subroutine which solves for X the overdetermined system AX = B of 
m linear equations in n unknowns for p right-hand sides. Written by Peter Businger, Computation 
Center, University of Texas (present address: Bell Telephone Laboratories, Murray Hill, N.J.), 
using the Businger-Golub algorithm [8|. 

MATH-PACK, ORTHLS, Orthogonal Polynomial Least-Squares Curve Fitting. One of the 
Univac 1108 MATH-PACK programs, written in FORTRAN V [40]. 

MPR3, Stepwise Multiple Regression with Variable Transformations. A FORTRAN II program 
written by 1Y1. A. Efroymson, Esso Research and Engineering Co., Madison, N.J., using the Efroym- 
son algorithm [16]. Available in the SHARE library: 7090- (2 3145MPR3 [26]. 

OMNITAB, a general-purpose computer program for statistical and numerical analysis. Devel- 
oped at the National Bureau of Standards by Joseph Hilsenrath et al. [23]. Now available in an 
A. S. A. FORTRAN version, OMNITAB allows the user to communicate with a computer in an 
efficient manner by means of simple English sentences. 

ORTHO. A program written by Philip J. Walsh, National Bureau of Standards (present address: 
University Computing Co., East Brunswick, N.J.), which uses a Gram-Schmidt orthonormaliza- 
tion process for least squares curve fitting. ORTHO exists as an ALGOL procedure [42], a FOR- 
TRAN program, a BASIC program (see LSFITW*** above), and as a routine of OMNITAB [23], 
where it is called by the commands FIT and POLYFIT. 

ORTHOL. A modification of the Davis-Rabinowitz orthonormalization algorithm [12, 13, 14], 
written in FORTRAN II by James W. Longley, Bureau of Labor Statistics, Washington, D.C., 
and Roger A. Blau, Bureau of Labor Statistics and Carnegie-Mellon University, Pittsburgh, Pa. 
[30]. 

POLFIT. An anonymous program written in BASIC for least squares polynomial curve fitting 
using orthogonal polynomials. 

POLRG, Polynomial Regression. One of the programs of the IBM System/360 Scientific Subroutine 
Package written in FORTRAN IV [24, 25]. 

SIMEX-***. A program written in BASIC for solving n simultaneous equations in n unknowns. 
Origin: Naval Ordnance Laboratory, Silver Spring, Md. Available in the C-E-I-R Mult i- Access 
Computer Services library [10]. 

SOLVER. A FORTRAN program written by Morris Newman, National Bureau of Standards, for 
obtaining the exact solution of the system 4X = B* or the inverse of a matrix A, by congruential 
methods [35]. The elements of A and B must be integers. 
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SPYMTX. A single precision FORTRAN IV program for inverting a matrix or solving a set of 
linear equations. To a program from the SHARE library (7090-F1 3180INV1 Single Precision 
Matrix Inversion with Selective Pivoting, written by A. R. Sadaka [26]), Sally T. Peavy, National 
Bureau of Standards, incorporated accuracy checks. 

STAT-PACK, GLH, General Linear Hypotheses. One of the Univac 1108 STAT-PACK programs, 
written in FORTRAN V [41]. 

STAT-PACK, REBSOM, Back Solution Multiple Regression. One of the Univac 1108 STAT- 
PACK programs, written in FORTRAN V [41]. 

STAT-PACK, RESTEM, Stepwise Multiple Regression. One of the Univac 1108 STAT-PACK 
programs, written in FORTRAN V [41]. 

STAT20***. A program written in BASIC for stepwise multiple linear regression. Written by 
Thomas E. Kurtz, Dartmouth College, Hanover, N.H. Available in the C-E-I-R Multi-Access Com- 
puter Services library [10]. 

STAT21***. A program written in BASIC for multiple linear regression with detailed output. 
Written by Gerald Childs, Dartmouth College, Hanover, N.H. Available in the C-E-I-R Multi-Access 
Computer Services library [10]. 

WRAP, Weighted Regression Analysis Program. A FORTRAN II program written by M. D. 
Fimple, Sandia Corp., Albuquerque, New Mexico. Available in the SHARE library: 7090-G2 3231 
WRAP [26]. 
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APPENDIX B 

DETAILS FOR TABLE 1 — SINGLE PRECISION (8 DIGITS) 



ALSQ 

BETA-HAT (Yi; 



COUNT 



1.0228963 


1.640 


.99553592 


2.350 


.99941321 


3.232 


1.0001108 


3.955 


.99999613 


5.412 


1.0000000 


8.000 




AVERAGE = 4.098 


BMD02R 




BETA-HAT (Yl) 


COUNT 


-17.13281 


-1.258 


39 . 34436 


-1.584 


-13.26675 


-1 .154 


2.92344 


-.284 


.89241 


.968 


1.00212 


2.674 




AVERAGE = -.106 


BMD03R 




BETA-HAT (Yl) 


COUNT 


394.23438 


-1 .898 


-16.00000 


-1 .230 


28.00000 


-1 .431 


-1 .00000 


-.301 


1 .00000 


6.000 


1.00049 


3.310 




AVERAGE = .742 


BMD03R 




BETA-HAT (Yl) 


COUNT 


5161 .95310 


-2.642 


40.000000 


-1 .591 


4.0000001 


-.477 


.50000000 


.301 


.89062500 


.961 


.99804688 


2.709 





1108 


EXAMPLE 1 


BETA-HAT (Y2) 




COUNT 


1.0000002 




6.699 


.10000037 




5.432 


.0099998263 




4.760 


.0010000220 




4.658 


.000099998902 




4.959 


.000010000020 




5.699 




AVERAGE = 


5.368 




1108 


EXAM r LE 2 


BETA-HAT (Y2) 




COUNT 


.99954 




3.337 


.10098775 




2.C05 


.0096306379 




1.433 


.0010499504 




1.301 


.000097199970 




1.553 


.000010055379 




2.257 




AVERAGE = 


1.981 




7094 


EXAMPLE 3 


BETA-HAT (Y2) 




COUNT 


1.04353 




1.361 


.10083 




2.081 


.01013 




1.886 


.00101 




2.000 


.00010 




2.000 


.00001 




1.000 




AVERAGE = 


1.721 




1108 


EXAMPLE 4 


BETA-HAT (Y2) 




COUNT 


1.07353 




1.134 


.10131836 




1.880 


.010009766 




3.010 


.00099945068 




3.260 


.000097751617 




1.648 


.0000099837780 




2.790 





AVERAGE = 


-.123 




AVERAGE = 


DAM 
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BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




2.20 




-.079 


1 .000 




.460 




.268 


.101 




.920 




1 .097 


.00975 




1.03 




1.523 


.00103 




.997 




2.523 


.0000982 




1.00 




3.000 


.0000100 






AVERAGE = 


1 .389 




AVERAGE = 


DAM 








1108 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




26.798895 




-1 .412 


1 .0000993 




-53.926606 




-1 .740 


.099779484 




21 .511053 




-1 .312 


.010084331 




-1 .7723664 




-.443 


.00098840467 




1 .1553755 




.809 


.00010065825 




.99692726 




2.512 


.0000099868524 






AVERAGE = 


-.264 




AVERAGE = 
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2.287 

EXAMPLE 

COUNT 
4.000 
2.000 
1.602 
1.523 
1.745 
3.000 

2.312 

EXAMPLE 

COUNT 

4.003 

2.657 

2.074 

1.936 

2.182 

2.881 

2.622 



DETAILS FOR TABLE 1 — SINGLE PRECISION (8 DIGITS) 



LINFIT (MILLER) 
BETA-HAT (Yl) 

7360.000 

-16598.000 

6379.500 

-877.906 

50.989 

.000 



7094 



COUNT 


BETA-HAT 


3.867 


1 .074 


4.220 


-.066 


3.805 


.074 


2.944 


-.008 


1 699 


.001 


.000 


.000 



(Y2) 



EXAMPLE 7 
COUNT 
1.131 
-.220 
-.806 
-.954 
-.954 
.000 





AVERAGE 


= -2.756 






AVERAGE = 


-.301 


LSTSQ 
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EXAMPLE 8 


BETA-HAT (Yl) 




COUNT 




BETA-HAT (Y2) 




COUNT 


.99973875 




3.583 




.99999997 




7.523 


1.0006891 




3.162 




.10000011 




5.959 


.99970413 




3.529 




.0099999484 




5.287 


1 .0000452 




4.345 




.0010000083 




5.081 


.99999718 




5.550 




.000099999460 




5.268 


1.0000001 




7.000 




.000010000012 




5.921 




AVERAGE 


= 4.528 






AVERAGE - 


5.840 


MATH-PACK 13.5, ORTHLS , ORTHOGONAL POLYNOMIAL CURVE FITTING 1108 


EXAMPLE 9 


BETA-HAT (Yl) 




COUNT 




BETA-HAT ( Y2 ) 




COUNT 


.94458008 




1.256 




.99999909 




6.041 


1 .1799316 




.745 




.10000322 




4.492 


.91607666 




1.076 




.0099984696 




3.815 


1 .0135651 




1.868 




.0010002495 




3.603 


.99912310 




3.057 




.000099983743 




3.789 


1.0000196 




4.708 




.000010000364 




4.439 




AVERAGE 


= 2.118 






AVERAGE = 


4.363 


MPR3, STEPWISE MULTIPLE REGRESSION, 


SHARE 


LIBRARY 3145MPR3 


7094 


EXAMPLE 10 


BETA-HAT (Yl) 




COUNT 




BETA-HAT (Y2) 




COUNT 


-20.24219 




-1.327 




.99933 




3.174 


43.49164 




-1.628 




.1013436 




1 .872 


-14.37052 




-1.187 




.009510736 




1 .310 


3.03207 




-.308 




.001065033 




1.187 


.88800 




.951 




.00009639953 




1 .444 


1.00219 




2.660 




.00001007054 




2.152 




AVERAGE 


= -.140 






AVERAGE - 


1 .856 


OMNITAB, USING THE 


MATRIX COMMANDS MTRANS , 


MMULT, INVERT 


7094 


EXAMPLE 11 


BETA-HAT (Yl) 




COUNT 




BETA-HAT (Y2) 




COUNT 


-91.999999 




-1.968 




.99780273 




2.658 


136.00000 




-2.130 




.10278320 




1 .555 


-48.000000 




-1.690 




.0086669922 




.875 


5.5000000 




-.653 




.0011596680 




.797 


.75000000 




.602 




.000092506409 




1 .125 


1.0063477 




2.197 




.000010177493 




1 .751 




AVERAGE 


= -.607 






AVERAGE = 


1 .460 


OMNITAB, USING THE 


MATRIX COMMANDS MTRANS, 


MMULT, INVERT 


1108 


EXAMPLE 12 


BETA-HAT (Yl) 




COUNT 




BETA-HAT (Y2) 




COUNT 


-220.65660 




-2.346 




.99512554 




2.312 


316.45422 




-2.499 




.10715805 




1 .145 


-96.578077 




-1.989 




.0077479167 




.647 


10.859373 




-.994 




.0012353163 




.628 


.52858572 




.327 




.000088502186 




.939 


1.0087148 




2.060 




.000010213975 




1 .670 




AVERAGE 


= -.907 






AVERAGE = 


1 .224 








82 









DETAILS FOR TABLE 1 



SINGLE PRECISION (8 DIGITS) 



OMNITAB. USING ORTHO SUBROUTINE 




7094 


EXAMPLE 


13 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 




1 .0012817 




2.892 


.99999949 




6.292 




.99780273 




2.658 


.10000013 




5.886 




.99932861 




3.173 


.0099999756 




5.613 




1 .0001755 




3.756 


.0010000018 




5.745 




.99998569 




4.844 


.000099999866 




5.873 




1 .0000004 




6.398 


.000010000004 




6.398 






AVERAGE = 


3.954 




AVERAGE = 


5.968 




OMNITAB, USING ORTHO SUBROUTINE 




1108 


EXAMPLE 


14 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 




1.0064697 




2.189 


.99999990 




7.000 




.99902344 




3.010 


.099999700 




5.523 




.99975586 




3.612 


.010000125 




4.903 




.99996948 




4.515 


.00099998200 




4.745 




1.0000100 




5.000 


.00010000109 




4.963 




.99999968 




6.495 


.0000099999778 




5.654 






AVERAGE = 


4.137 




AVERAGE = 


5.464 




ORTHO, WITH RE-ORTHOGONALIZATION OMITTED 




1108 


EXAMPLE 


15 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 




-1216.5426 




-3.085 


.98419483 




1.801 




2752.0557 




-3.439 


.13523918 




.453 




-1057.0931 




-3.025 


-.0034660707 




-.129 




146.97336 




-2.164 


.0028495983 




-.267 




-7.3080225 




-.919 


-.0000049256487 




-.021 




1 . 1663037 




.779 


.000012094996 




.679 






AVERAGE = 


-1 .976 




AVERAGE = 


.419 




ORTHOL 








1108 


EXAMPLE 


16 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 




. 99784447 




2.666 


1.0000000 




8.000 




.98687472 




1 .882 


.099999778 




5.654 




1.0029743 




2.527 


.010000041 




5.387 




.99961372 




3.413 


.00099999654 




5.461 




1.0000213 




4.672 


.00010000013 




5.886 




.99999960 




6.398 


.0000099999984 




6.796 






AVERAGE = 


3.593 




AVERAGE = 


6.197 




POLRG, IBM SYSTEM/360 SCIENT 


IFIC SUBROUTINE PACKAGE 


1108 


EXAMPLE 


17 


BETA-HAT (Yl) 




COUNT 


BETA-HAT ( Y2 ) 




COUNT 




-1823.8047 




-3.261 


.98438931 




1.807 




28.622013 




-1.441 


.10009000 




3.046 




-3.6844511 




-.671 


.010146444 




1.834 




3.0450442 




-.311 


.0010054175 




2.266 




.93157484 




1.165 


.000099141363 




2.066 




1.0004238 




3.373 


.000010021910 




2.659 






AVERAGE = 


-.191 




AVERAGE = 


2.280 




SPVMTX 








1108 


EXAMPLE 


18 


BETA-HAT (Yl) 




COUNT 


BETA-HAT ( Y2 ) 




COUNT 




64.191025 




-1.801 


1.0014403 




2.842 




-134.65426 




-2.132 


.097101737 




1.538 




51.859977 




-1.706 


.011047948 




.980 




-5.8942945 




-.838 


.00086162069 




.859 




1.3872223 




.412 


.00010761766 




1.118 




.99232943 




2.115 


.0000098514820 




1.828 






AVERAGE = 


-.658 




AVERAGE = 


1 .527 
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DETAILS FOR TABLE 1 



SINGLE PRECISION (8 DIGITS) 



STAT-PACK 8.13, GLH , GENERAL LINEAR HYPOTHESIS 



1108 



BETA-HAT (Yl) 




COUNT 


-11.970093 




-1.113 


27.245361 




-1.419 


-8.5661011 




-.981 


2.2717514 




-.104 


.92961073 




1.152 


1.0013784 




2.861 




AVERAGE = 


.066 



BETA-HAT (Y2) 
.99989064 
.10018034 
.0099404901 
.0010073970 
.000099610348 
.000010007344 



AVERAGE 



EXAMPLE 19 

COUNT 

3.961 

2.744 

2.225 

2.131 

2.409 

3.134 

2.767 



STAT-PACK 9.2, REBSOM, BACK SOLUTION MULTIPLE REGRESSION 

BETA-HAT (Yl) COUNT BETA-HAT (Y2) 

-11.97009 -1.113 .9998906 

27.24536 -1.419 .1001803 

-8.56610 -.981 .009940490 

2.27175 -.104 .001007397 

.92961 1.152 .00009961035 

1.00138 2.860 .00001000734 



1108 



EXAMPLE 20 

COUNT 

3.961 

2.744 

2.225 

2.131 

2.409 

3.134 





AVERAGE 


.066 


T-PACK 9.1, RESTEM, 


STEPWISE MUL' 


BETA-HAT (Yl) 






COUNT 


-1.5703125 






-.410 


7.5583214 






-,817 


-1.5695286 






-.410 


1.3551083 






.450 


.97985181 






1.696 


1.0004015 






3.396 




AVERAGE 


.651 


lP, WEIGHTED REGRESSION, 


SHARE LI 


BETA-HAT (Yl) 






COUNT 


2991622. 






-6 . 476 


-6065892. 






-6.783 


2218821. 






-6.346 


-296194.5 






-5.472 


16462.20 






-4.216 


-322.5731 






-2.510 




AVERAGE 


= -5.300 



AVERAGE 



1108 



BETA-HAT (Y2) 

1.0002102 
.099611598 
.010136487 
.00098222795 
.00010097076 
.0000099811599 



BETA-HAT (Y2) 
-33.84546 
71.54880 
-26.16913 
3.493256 
-.1936966 
.003812985 



7094 



2.767 

EXAMPLE 21 

COUNT 

3.677 

.411 

.865 

.750 

.013 



2. 

1. 
1. 
2. 
2.725 



AVERAGE = 2.407 



EXAMPLE 22 

COUNT 
-1.542 
-2.854 
-3.418 
-3.543 
-3.287 
-2.580 



AVERAGE 



-2.871 



DETAILS FOR TABLE 2 — SINGLE PRECISION (9 DIGITS) 



LINFIT 



BETA-HAT (Yl) 




COUNT 


2.76639 




-.247 


-2.72473 




-.571 


2 . 38633 




-.142 


.812925 




.728 


1.01048 




1.980 


.999793 




3.684 




AVERAGE = 


.905 


F 

BETA-HAT (Yl) 




COUNT 


-5.5 




-.813 


-12. 




-1.114 


-3. 




-.602 


0. 




.000 


.957031 




1.367 


.999023 




3.010 




AVERAGE = 


.308 





235 


EXAMPLE 1 


BETA-HAT (Y2) 




COUNT 


1.00006 




4.222 


.0998764 




2.908 


.010045 




2.347 


.000994014 




2.223 


.000100332 




2.479 


.00000999350 




3.187 




AVERAGE = 


2.894 




235 


EXAMPLE 2 


BETA-HAT (Y2) 




COUNT 


.999893 




3.971 


.099762 




2.623 


.00991058 




2.049 


.000970840 




1.535 


.0000993013 




2.156 


.00000997260 




2.562 




AVERAGE = 


2.483 



84 



DETAILS FOR TABLE 2 



SINGLE PRECISION (9 DIGITS) 



LSFITW 










235 


EXAMPLE 


BETA HAT (Yl) 




COUNT 


BETA-HAT (Y2) 






COUNT 


.999130249 




3.061 


.999999897555 






6.990 


.99761963 




2.623 


.0999999824213 






6.755 


1.00102997 




2.987 


.0100000116561 






5.933 


.999854088 




3.836 


.00099999815736 




5.735 


1.00000715256 




5.146 


.000100000104819 




5.980 


.999999890104 




6.959 


.00000999999813838 


6.730 




AVERAGE = 


4.102 




AVERAGE = 


6.354 


POLFIT 










235 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 






COUNT 


.99387360 




2.213 


.9999999618158 






7.418 


1.00894165 




2.049 


.1000000573928 






6.241 


.99534607 




2.332 


.0099999622960 






5.424 


1.000703812 




3.153 


.00100000655382 




5.184 


.9999548197 




4.345 


.000099999526381 




5.325 


1.000000998378 




6.001 


.0000100000114824 




5.940 




AVERAGE = 


3.349 




AVERAGE = 


5 .922 


SIMEX 










235 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 






COUNT 


1.74226 




.129 


.999966 






4.469 


-.313568 




-.118 


.100063 






3.201 


1 .44267 




.354 


.00997838 






2.665 


.944463 




1 .255 


.00100276 






2.559 


1.00294 




2.532 


.0000998520 






2.830 


.999945 




4.260 


.0000100028 






3.553 




AVERAGE = 


1 402 




AVERAGE = 


3.213 


STAT20 










235 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 






COUNT 


4.70801 




-.569 


1.00006 






4.222 


-6.48121 




-.874 


.0998837 






2.934 


3.72065 




-.435 


.010042 






2.377 


.638874 




.442 


.000994436 






2.255 


1 .01997 




1 .700 


.000100307 






2.513 


.999609 




3.408 


.00000999400 






3.222 




AVERAGE = 


.612 




AVERAGE = 


2.920 


STAT21 

BETA-HAT (Yl) 

2.089 
-1.11166 
1 .75217 

.901511 
1.00539 
.999895 




COUNT 
-.037 
-.325 
.124 
1.007 
2.268 
3.979 


BETA-HAT (Y2) 

1.00003 
.0999349 
.0100234 
.000996913 
.000100170 
00000999668 




235 


EXAMPLE 

COUNT 

4.523 

3.186 

2.631 

2.510 

2.770 

3.479 




AVERAGE = 


1 .169 




AVERAGE = 


3.183 
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DETAILS FOR TABLE 3 — DOUBLE PRECISION (16 DIGITS) 



BMD05R 



BETA-HAT (Yl) 


COUNT 


1.0000003 


6.523 


.99999920 


6.097 


1.0000002 


6.699 


.99999996 


7.398 


.9999999 


7.000 


.99999999 


8.000 


AVERAGE = 


6.953 


DPVMTX 




BETA-HAT (Yl) 


COUNT 


1.000000206882520 


6.684 


.9999995965200948 


6.394 


1.000000145134836 


6.838 


.9999999808270597 


7.717 


1.000000001057576 


8.976 


.9999999999793303 


10.685 


AVERAGE = 


7.882 



7094 


EXAMPLE 1 


BETA-HAT (Y2) 


COUNT 


.99999998 


7.699 


.10000004 


6.398 


.0099999798 


5.695 


.0010000031 


5.509 


.000099999792 


5.682 


.000010000004 


6.398 


AVERAGE = 


6.230 


1107 


EXAMPLE 2 


BETA-HAT (Y2) 


COUNT 


1.000000000006017 


11.221 


.09999999998899705 


9.958 


.01000000000383598 


9.416 


. 0009999999995039349 


9.304 


. 0001000000000269390 


9.570 


. 000009999999999479781 


10.284 


AVERAGE = 


9.959 



DETAILS FOR TABLE 4 — DOUBLE PRECISION (18 DIGITS) 



BETA-HAT (Yl) 


COUNT 


1.000000000005630 


11.249 


.9999999999933983 


11.180 


1.000000000002151 


11.667 


.9999999999997248 


12.560 


1.000000000000015 


13.824 


.9999999999999997 


15.520 


AVERAGE 


= 12.667 


BMD02R 




BETA-HAT (Yl) 


COUNT 


. 9999999968749762 


8.505 


1.000000006749235 


8.171 


.9999999974683392 


8.597 


1 .000000000342994 


9.465 


.9999999999807494 


10.716 


1 .000000000000381 


12.419 


AVERAGE 


= 9 . 645 


BMD05R 




BETA-HAT (Yl) 


COUNT 


1.00000000634827302 


8.197 


. 999999987016701379 


7.887 


1.00000000476533960 


8.322 


. 999999999362724245 


9.196 


1 .00000000003545420 


10.450 


.999999999999302617 


12.157 


AVERAGE 


= 9 . 368 


DPVMTX 




BETA-HAT (Yl) 


COUNT 


.9999999971390853 


8.543 


1.000000005557233 


8.255 


.9999999980049357 


8.700 


1.000000000263132 


9.580 


.9999999999855033 


10.839 


1.000000000000283 


12.548 


AVERAGE 


= 9 . 744 



1108 

BETA-HAT (Y2) 

1.00000000000000000 
. 0999999999999999140 
.0100000000000000339 



EXAMPLE 1 

COUNT 
18.000 
15.065 
14.470 



. 000999999999999995720 14 . 368 



. 000100000000000000224 



14.650 



.00000999999999999999580 15.376 



AVERAGE = 15.322 



1108 


EXAMPLE 


BETA-HAT (Y2) 


COUNT 


1 .000000000000007 


14.155 


.09999999999998616 


12.859 


.01000000000000481 


12.318 


. 0009999999999993756 


12.205 


.0001000000000000339 


12.470 



. 000009999999999999342 13 . 182 



AVERAGE = 12.865 



1108 


EXAMPLE 3 


BETA-HAT (Y2) 


COUNT 


1 .00000000000007691 


13.114 


. 0999999999998442950 


11.808 


.0100000000000568478 


11 .245 


. 000999999999992425020 


11 .121 


.000100000000000420311 


11 .376 


. 00000999999999999174930 


12.084 


AVERAGE = 


11.791 


1108 


EXAMPLE 4 


BETA-HAT (Y2) 


COUNT 


.9999999999999998 


15.693 


.1000000000000033 


13.481 


.0099999999999983 r 


12.790 


.0010000000000002' 


12.607 


.00009999999999998 7 


12.826 


.00001000000000000 1 


13.509 


AVERAGE = 


13.484 
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DETAILS FOR TABLE 4 



DOUBLE PRECISION (18 DIGITS) 



LSTSQ 




1108 


EXAMPLE 5 


BETA-HAT (Yl) 


COUNT 


BETA-HAT (Y2) 


COUNT 


. 999999999999999464 


15.269 


1.00000000000000000 


18.000 


. 999999999999943084 


13.245 


. 0999999999999999950 


16.284 


1 .00000000000004282 


13.368 


. 0100000000000000021 


15.683 


. 999999999999991067 


14.049 


. 000999999999999999710 


15.538 


1 .00000000000000068 


15.169 


. 000100000000000000017 


15.771 


. 999999999999999985 


16.761 


. 00000999999999999999970 


16.480 


AVERAGE 


= 14.643 


AVERAGE = 


16.293 


MATH-PACK 13.5, ORTHLS , ORTHOGONAL POLYNOMIAL CURVE FITTING 1108 


EXAMPLE 6 


BETA-HAT (Yl) 


COUNT 


BETA-HAT (Y2) 


COUNT 


. 999999999992709829 


11.137 


.999999999999999901 


15.997 


1 .00000000001946887 


10.711 


.100000000000000225 


14.648 


. 999999999991466379 


11 .069 


. 00999999999999988980 


13.958 


1 .00000000000133404 


11.875 


.00100000000000001872 


13.728 


. 999999999999915178 


13.071 


. 0000999999999999987360 


13.898 


1 .00000000000000188 


14.727 


. 0000100000000000000292 


14.535 


AVERAGE 


= 12.098 


AVERAGE = 


14.461 


ORTHO 




1108 


EXAMPLE 7 


BETA-HAT (Yl) 


COUNT 


BETA-HAT (Y2) 


COUNT 


.999999999997051246 


11 .530 


. 999999999999999946 


16.242 


1 .00000000000051159 


12.291 


. 0999999999999999970 


16.488 


1 .00000000000034817 


12.458 


.0100000000000000078 


15. 109 


. 999999999999897859 


12.991 


. 000999999999999998330 


14.777 


1 .00000000000000760 


14.119 


.000100000000000000121 


14.918 


. 999999999999999820 


15.740 


. 00000999999999999999720 


15.550 


AVERAGE 


= 13.188 


AVERAGE = 


15.514 


ORTHO. WITH RE-ORTHOGONALIZATION OMITTED 


1108 


EXAMPLE 8 


BETA-HAT (Yl) 


COUNT 


BETA-HAT (Y2) 


COUNT 


.999999858679196051 


6.850 


.999999999998151750 


11.733 


1 .00000031785242527 


6.498 


.100000000004101749 


10.387 


. 999999878054865120 


6.914 


. 00999999999843660520 


9.806 


1 .00000001679285067 


7.775 


. 00100000000021433129 


9.669 


.999999999045590670 


9.020 


. 0000999999999878591650 


9.916 


1 .00000000001908297 


10.719 


. 0000100000000002421209 


10.616 



AVERAGE = 7.963 



AVERAGE = 10.354 



COUNT 


12 


610 


11 


592 


12 


086 


12 


948 


14 


175 


15 


863 



ORTHOL 

BETA-HAT (Yl) 

. 999999999999754768 
1 .00000000000255572 

. 999999999999179468 
1 .00000000000011276 

. 999999999999993318 
1 .00000000000000014 



AVERAGE = 13.212 

POLRG, IBM SYSTEM/360 SCIENTIFIC SUBROUTINE PACKAGE 

BETA-HAT (Yl) COUNT 

1 .00000011110114428 6.954 

. 999999990457514712 8 . 020 

1 . 00000000321449412 8 . 493 

. 999999999544922239 9 .-342 

1 . 00000000002322160 10 . 634 



.999999999999491835 



AVERAGE 



12.294 



.290 



1108 


EXAMPLE 9 


BETA-HAT (Y2) 


COUNT 


1.00000000000000000 


18.000 


.100000000000000036 


15.444 


. 00999999999999998660 


14.873 


. 00100000000000000191 


14.719 


. 0000999999999999998880 


14.950 


. 0000100000000000000023 


15.640 


AVERAGE = 


15.604 


PACKAGE 1108 


EXAMPLE 10 


BETA-HAT (Y2) 


COUNT 


1.00000000000136645 


11.864 


. 099999999999941 1990 


12.231 


. 0100000000000305940 


11.514 


. 000999999999995547370 


11.351 


. 000100000000000227181 


11.644 


. 00000999999999999416480 


12.234 


AVERAGE = 


11 .806 



87 



DETAILS FOR TABLE 4 



DOUBLE PRECISION (18 DIGITS) 



STAT-PACK 9.1, RESTEM, STEPWISE MULTIPLE 

BETA-HAT (Yl) COUNT 

1 . 00000000453928806 8 . 343 

. 999999990379828084 8 . 017 

1 . 00000000358299878 8 . 446 

. 999999999516629162 9.316 

1 . 00000000002705048 10 . 568 

. 999999999999465664 12 . 272 



REGRESSION 1108 EXAMPLE 11 

BETA-HAT (Y2) COUNT 

1 . 00000000000004872 13 . 312 

. 0999999999999056810 12 . 025 

. 0100000000000336338 11 . 473 

. 000999999999995591090 11 . 356 

. 000100000000000241649 11 . 617 

.00000999999999999530180 12.328 



AVERAGE = 9.494 



AVERAGE = 12.019 



DETAILS FOR TABLE 5 



SINGLE PRECISION (8 DIGITS), WITH INNER PRODUCTS 
ACCUMULATED IN DOUBLE PRECISION (18 DIGITS) 



ALSQ 








1108 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 


1.0366969 




1.435 


1.0000001 




7.000 


.99202651 




2.098 


.099999869 




5.883 


.99865498 




2.871 


.010000017 




5.770 


1.0003119 




3.506 


.00099999901 




6.004 


.99998119 




4.726 


.00010000003 




6.523 


1.0000004 




6.398 


.0000099999999 




8.000 




AVERAGE = 


3.506 




AVERAGE = 


6.530 


LSTSQ 








1108 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 


1.0000000 




8.000 


.99999999 




8.000 


.99999999 




8.000 


.10000004 




6.398 


1.0000000 




8.000 


.0099999798 




5.695 


1.0000000 




8.000 


.0010000032 




5.495 


1.0000000 




8.000 


.000099999794 




5.686 


1.0000000 




8.000 


.000010000004 




6.398 




AVERAGE = 


8.000 




AVERAGE = 


6.279 


ORTHO 








1108 


EXAMPLE 


BETA-HAT (Yl) 




COUNT 


BETA-HAT (Y2) 




COUNT 


1.0027425 




2.562 


1.0000003 




6.523 


. 99674778 




2.488 


.099999898 




5.991 


1.0013667 




2.864 


.010000011 




5.959 


.99983171 




3.774 


.00099999974 




6.585 


1.0000096 




5.018 


.000099999980 




6.699 


.99999981 




6.721 


.000010000001 




7.000 




AVERAGE = 


3.904 




AVERAGE = 


6.459 



DETAILS FOR TABLE 6 



MULTIPLE PRECISION INTEGER ARITHMETIC 



SOLVER 

BETA-HAT 



[Yl) 



1.00000000000000000 
1.00000000000000000 
1 .00000000000000000 
1 .00000000000000000 
1 .00000000000000000 
1 .00000000000000000 



COUNT 
18.000 
18.000 
18.000 
18.000 
18.000 
18.000 



AVERAGE = 18.000 



1108 


EXAMPLE 1 


BETA-HAT (Y2) 


COUNT 


.999999999999999995 


17.159 


. 0999999999999999995 


17.187 


. 00999999999999999996 


17.169 


. 000999999999999999996 


17.294 


. 0000999999999999999996 


17.276 


. 0000100000000000000000 


18.000 


AVERAGE = 


17.347 
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